% ISOM_01_model
%%%%%%%%%%%%% Parameters %%%%%%%%%%%%%%%%



T=100000;
cut=500;  % Initial condition dependence
N=1;

addpath('auxiliar')
 load seed
 rng(s);
 slength=30; %Total Length of the Sample (used+Cut)
finalcut=17; % Periods Of Market access before we use (cut)


default_zone=zeros(NBG,NB,NSS);


 for j=1:NSS %Today's Exogenous Shock
           for t=2:NB %Today's private debt index;                     
                for i=1:NBG %Today's public Debt shock
                      if Def_func(i,t,j)< Def_func(i,t-1,j)  
                          default_zone(i,t,j)=1;
                      end
                 end
           end
 end


% 1. Simulation
S=zeros(N,T+cut,N);
S_index=zeros(T+cut,N);

inv_dist=invariant(Prob);
cum_dist=cumsum(inv_dist);
load seed
rng(s);
init0=rand(N,1);
%FInd Invariant Distribution
inva = @(r) find(r<cum_dist,1,'first'); % find the 1st index s.t. r<D(i);
% Now this are your results of the random trials
rR = arrayfun(inva,init0);

for i=1:N
rng(i,'twister');
S=markov(Prob,T+cut+1,rR(i),1:length(Prob));
S_index(:,i)=S';
end


ySIM=zeros(T+cut,N);
yNSIM=zeros(T+cut,N);
KAPPASIM=zeros(T+cut,N);
epsilons_default=zeros(T+cut,N);
epsilons_debt=zeros(T+cut,N);

rng(N+1,'twister');
epsilons_default=rand(T+cut,N);
epsilon_reentry=rand(T+cut,N);
rng(N+2,'twister');
epsilons_debt=rand(T+cut,N);
PVSIM=zeros(T+cut,N);
SSSIM=zeros(T+cut,N);
for i=1:T+cut
   for j=1:N
    ySIM(i,j)=yT(S_index(i,j));
    yNSIM(i,j)=yN(S_index(i,j));
    KAPPASIM(i,j)=KAPPAS(S_index(i,j));
    PVSIM(i,j)=PVD(S_index(i,j));
   end
end


% Core Interpolations

bpSIM=zeros(T+cut,N);
bSIM=zeros(T+cut,N);
BGSIM=zeros(T+cut,N);
BGpSIM=zeros(T+cut,N);
QSIM=zeros(T+cut,N);
DefSIM=zeros(T+cut,N);
BGp_bigSIM=zeros(T+cut,N,NBG);
qSIM=zeros(T+cut,N);
DebtLimSIM=zeros(T+cut,N);
BGSIM_Defaulted=zeros(T+cut,N);
bSIM_Defaulted=zeros(T+cut,N);
      TransSIM=zeros(T+cut,N);
            DistSIM=zeros(T+cut,N);
            DistSIM_DEFAULTED_4Q=zeros(T+cut,N);
            DistSIM_DEFAULTED_12Q=zeros(T+cut,N);
            Transfer_DEFAULTED_4Q=zeros(T+cut,N);
            Transfer_DEFAULTED_12Q=zeros(T+cut,N);

bSIM_Defaultedm1=zeros(T+cut,N);
bSIM_Defaultedm2=zeros(T+cut,N);
bSIM_Defaultedm3=zeros(T+cut,N);
bSIM_Defaultedm4=zeros(T+cut,N);
bSIM_Defaultedm5=zeros(T+cut,N);

BGSIM_Defaultedm1=zeros(T+cut,N);
BGSIM_Defaultedm2=zeros(T+cut,N);
BGSIM_Defaultedm3=zeros(T+cut,N);
BGSIM_Defaultedm4=zeros(T+cut,N);
BGSIM_Defaultedm5=zeros(T+cut,N);
reentrySIM=zeros(T+cut,N);

zoneSIM=NaN(T+cut,N);
DefSP_SIM=NaN(T+cut,N);


%Indices
BGiSIM=NaN(T+cut+1,N);
BGpiSIM=NaN(T+cut,N);
cSIM=zeros(T+cut,N);
priceSIM=zeros(T+cut,N);
SpreadSIM=zeros(T+cut,N);
zeroBGI_ind=find(BGgrid>-0.0001,1);
BadStandindSIM=zeros(T+cut,N);
BGiSIM(1,:)=zeroBGI_ind;
bSIM(1,:)=0;%0.5;

for i=1:T+cut
    for j=1:N
    bad_stand=BadStandindSIM(i,j);
    yindex=S_index(i,j);
    BGindex=BGiSIM(i,j);
    qSIM(i,j)=q(S_index(i,j));
    bgrid_cc_temp=squeeze(bgrid(BGindex,:));
    zoneSIM(i,j)=interp1(bgrid_def_temp,squeeze(Def_func(BGindex,:,yindex)),bSIM(i,j),'linear','extrap');
    def_temp=min(max(interp1(bgrid_cc_temp,squeeze(Def_func(BGindex,:,yindex)),bSIM(i,j),'linear','extrap'),0),1);
   if bad_stand==0 || (bad_stand==1 && theta>epsilon_reentry(i,j)) % WE ARE IN GOOD STANDING or we were in bad standing but we reenter
       if (bad_stand==1 && theta>epsilon_reentry(i,j)) % THIS WAS A REENTRY PERIOD
           reentrySIM(i,j)=1;
       end
        if def_temp>epsilons_default(i,j) %WE ENTER A DEFAULT
             def_temp_SP=min(max(interp1(bgrid_cc_temp,squeeze(Def_funcSP(BGindex,:,yindex)),bSIM(i,j),'linear','extrap'),0),1);
            if def_temp_SP>epsilons_default(i,j)
                DefSP_SIM(i,j)=0; % THE PLANNER WOULD HAVE ALSO ENTERED THIS DEFAULT
            else
                 DefSP_SIM(i,j)=1; % THE PLANNER WOULD NOT HAVE ALSO ENTERED THIS DEFAULT
            end

            BadStandindSIM(i+1,j)=1; % TOMORROW WILL BE IN BAD STANDING
            DefSIM(i,j)=true(1);            
            BGpSIM(i,j)=0;
            BGSIM(i,j)=0;
            BGSIM_Defaulted(i,j)=BGgrid(BGindex);
            bSIM_Defaulted(i,j)=bSIM(i,j);
            BGpiSIM(i,j)=zeroBG_ind;
            BGiSIM(i,j)=zeroBG_ind;
            BGindex=zeroBG_ind;
            bgrid_def_temp=squeeze(bgrid(BGindex,:));
            bpSIM(i,j)=interp1(bgrid_def_temp,squeeze(bp_bad(:,yindex)),bSIM(i,j),'linear','extrap');
            cSIM(i,j)=max(interp1(bgrid_def_temp,squeeze(c_big_bad(:,yindex)),bSIM(i,j),'linear','extrap'),1e-8);
           
            priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
            QSIM(i,j)=interp1(bgrid_def_temp,squeeze(Qgov(zeroBG_ind,:,yindex,zeroBG_ind)),bSIM(i,j),'linear','extrap');
            SpreadSIM(i,j)=0;
            DebtLimSIM(i,j)=(KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))/qSIM(i,j);
            bpSIM(i,j)=min(bpSIM(i,j), DebtLimSIM(i,j));
            TransSIM(i,j)=QSIM(i,j)*(BGpSIM(i,j)-(1-delta)*BGSIM(i,j))-delta*BGSIM(i,j);
        


            DistSIM(i,j)=-(bpSIM(i,j)-(1/qSIM(i,j))*KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)));        
            DistSIM_DEFAULTED_4Q(i,j)=mean(DistSIM((i-5):i-1,j));
            DistSIM_DEFAULTED_12Q(i,j)=mean(DistSIM((i-1):i-1,j));
            Transfer_DEFAULTED_4Q(i,j)=mean(TransSIM((i-5):i-1,j));
            Transfer_DEFAULTED_12Q(i,j)=mean(TransSIM((i-1):i-1,j));

            bSIM_Defaultedm1(i,j)=bSIM(i-1,j);
            bSIM_Defaultedm2(i,j)=bSIM(i-2,j);
            bSIM_Defaultedm3(i,j)=bSIM(i-3,j);
            bSIM_Defaultedm4(i,j)=bSIM(i-4,j);
            bSIM_Defaultedm5(i,j)=bSIM(i-5,j);
            BGSIM_Defaultedm1(i,j)=BGSIM(i-1,j);
            BGSIM_Defaultedm2(i,j)=BGSIM(i-2,j);
            BGSIM_Defaultedm3(i,j)=BGSIM(i-3,j);
            BGSIM_Defaultedm4(i,j)=BGSIM(i-4,j);
            BGSIM_Defaultedm5(i,j)=BGSIM(i-5,j);


        else
            DefSIM(i,j)=false(1);  
            BadStandindSIM(i+1,j)=0; % TOMORROW WILL BE IN GOOD STANDING
            for jj=1:NBG
                gtemp=min(max(interp1(bgrid_cc_temp,squeeze(G_big(BGindex,:,yindex,jj)),bSIM(i,j),'linear','extrap'),0),1);
                BGp_bigSIM(i,j,jj)=gtemp+ BGp_bigSIM(i,j,max(jj-1,1)); %We build the Cumulative Distribution
                if epsilons_debt(i,j)<=BGp_bigSIM(i,j,jj) && isnan(BGpiSIM(i,j))==1 % We find the spot in which the shoc crosses the cumulative
                    BGpiSIM(i,j)=jj;
                end
            end             
            BGpSIM(i,j)=BGgrid(BGpiSIM(i,j));
            BGSIM(i,j)=BGgrid(BGindex);
            bpSIM(i,j)=interp1(bgrid_cc_temp,squeeze(bp(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap');
            cSIM(i,j)=max(interp1(bgrid_cc_temp,squeeze(c_big(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap'),1e-8);
            priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
           
            QSIM(i,j)=interp1(bgrid_cc_temp,squeeze(Qgov(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap');
            SpreadSIM(i,j)=(1+(delta)/QSIM(i,j)-delta)-(1+rbar);
            DebtLimSIM(i,j)=(KAPPASIM(i,j)*((priceSIM(i,j)*yNSIM(i,j))+ySIM(i,j)))/qSIM(i,j);
            bpSIM(i,j)=min(bpSIM(i,j), DebtLimSIM(i,j));
                        BGSIM_Defaulted(i,j)=NaN;
                        bSIM_Defaulted(i,j)=NaN;
             TransSIM(i,j)=QSIM(i,j)*(BGpSIM(i,j)-(1-delta)*BGSIM(i,j))-delta*BGSIM(i,j);
            DistSIM(i,j)=-(bpSIM(i,j)-(1/qSIM(i,j))*KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)));        
            DistSIM_DEFAULTED_4Q(i,j)=NaN;
            DistSIM_DEFAULTED_12Q(i,j)=NaN;
            Transfer_DEFAULTED_4Q(i,j)=NaN;
            Transfer_DEFAULTED_12Q(i,j)=NaN;
            bSIM_Defaultedm1(i,j)=NaN;
            bSIM_Defaultedm2(i,j)=NaN;
            bSIM_Defaultedm3(i,j)=NaN;
            bSIM_Defaultedm4(i,j)=NaN;
            bSIM_Defaultedm5(i,j)=NaN;
            BGSIM_Defaultedm1(i,j)=NaN;
            BGSIM_Defaultedm2(i,j)=NaN;
            BGSIM_Defaultedm3(i,j)=NaN;
             BGSIM_Defaultedm4(i,j)=NaN;
            BGSIM_Defaultedm5(i,j)=NaN;


        end

   else % WE ARE IN BAD STANDING and we dont enter
          BadStandindSIM(i+1,j)=1; % TOMORROW WILL BE IN BAD STANDING
            DefSIM(i,j)=true(1);            
            BGpSIM(i,j)=0;
            BGSIM(i,j)=0;
            BGSIM_Defaulted(i,j)=NaN; %BGgrid(BGindex);
            bSIM_Defaulted(i,j)=NaN;
            BGpiSIM(i,j)=zeroBG_ind;
            BGiSIM(i,j)=zeroBG_ind;
            BGindex=zeroBG_ind;
            bgrid_def_temp=squeeze(bgrid(BGindex,:));
            bpSIM(i,j)=interp1(bgrid_def_temp,squeeze(bp_bad(:,yindex)),bSIM(i,j),'linear','extrap');
            cSIM(i,j)=max(interp1(bgrid_def_temp,squeeze(c_big_bad(:,yindex)),bSIM(i,j),'linear','extrap'),1e-8);
            %priceSIM(i,j)=interp1(bgrid_def_temp,squeeze(price(zeroBG_ind,:,yindex,zeroBG_ind)),bSIM(i,j),'linear','extrap');
            priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
            QSIM(i,j)=interp1(bgrid_def_temp,squeeze(Qgov(zeroBG_ind,:,yindex,zeroBG_ind)),bSIM(i,j),'linear','extrap');
            SpreadSIM(i,j)=0;
            DebtLimSIM(i,j)=(KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))/qSIM(i,j);
            bpSIM(i,j)=min(bpSIM(i,j), DebtLimSIM(i,j));
            TransSIM(i,j)=QSIM(i,j)*(BGpSIM(i,j)-(1-delta)*BGSIM(i,j))-delta*BGSIM(i,j);
            DistSIM(i,j)=-(bpSIM(i,j)-(1/qSIM(i,j))*KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)));        
            DistSIM_DEFAULTED_4Q(i,j)=NaN;
            DistSIM_DEFAULTED_12Q(i,j)=NaN;
            Transfer_DEFAULTED_4Q(i,j)=NaN;
            Transfer_DEFAULTED_12Q(i,j)=NaN;
             bSIM_Defaultedm1(i,j)=NaN;
            bSIM_Defaultedm2(i,j)=NaN;
            bSIM_Defaultedm3(i,j)=NaN;
            bSIM_Defaultedm4(i,j)=NaN;
            bSIM_Defaultedm5(i,j)=NaN;
            BGSIM_Defaultedm1(i,j)=NaN;
            BGSIM_Defaultedm2(i,j)=NaN;
            BGSIM_Defaultedm3(i,j)=NaN;
            BGSIM_Defaultedm4(i,j)=NaN;
            BGSIM_Defaultedm5(i,j)=NaN;
            

   end
    BGiSIM(i+1,j)=BGpiSIM(i,j);
    bSIM(i+1,j)=bpSIM(i,j);
    
    end
end

% Computation

bSIM=bSIM(1:end-1,:);%

% Current Accpunt Calculation
CA_SIM=bpSIM(2:end,:)-bpSIM(1:end-1,:);  % Current Account of the Private Sector

cSIM=cSIM(cut+1:T+cut,:);
ySIM=ySIM(cut+1:T+cut,:);
yNSIM=yNSIM(cut+1:T+cut,:);
S_index=S_index(cut+1:T+cut,:);
DefSIM=DefSIM(cut+1:T+cut,:); 
BGiSIM=BGiSIM(cut+1:T+cut,:);
BGpiSIM=BGpiSIM(cut+1:T+cut,:); 
BGpSIM=BGpSIM(cut+1:T+cut,:); 
bpSIM=bpSIM(cut+1:T+cut,:);
BGSIM=BGSIM(cut+1:T+cut,:);
bSIM=bSIM(cut+1:T+cut,:);
QSIM=QSIM(cut+1:T+cut,:);
PVSIM=PVSIM(cut+1:T+cut,:);
KAPPASIM=KAPPASIM(cut+1:T+cut,:);
BGp_bigSIM=BGp_bigSIM(cut+1:T+cut,:,:);
epsilons_debt=epsilons_debt(cut+1:T+cut,:);
epsilons_default=epsilons_default(cut+1:T+cut,:);
DebtLimSIM=DebtLimSIM(cut+1:T+cut,:);
qSIM=qSIM(cut+1:T+cut,:);
reentrySIM=reentrySIM(cut+1:T+cut,:);
DefSP_SIM=DefSP_SIM(cut+1:T+cut,:);

BGSIM_Defaulted=BGSIM_Defaulted(cut+1:T+cut,:);
bSIM_Defaulted=bSIM_Defaulted(cut+1:T+cut,:);
     TransSIM=TransSIM(cut+1:T+cut,:);
     
            DistSIM_=DistSIM(cut+1:T+cut,:);
            DistSIM_DEFAULTED_4Q=DistSIM_DEFAULTED_4Q(cut+1:T+cut,:);
            DistSIM_DEFAULTED_12Q=DistSIM_DEFAULTED_12Q(cut+1:T+cut,:);
            Transfer_DEFAULTED_4Q=Transfer_DEFAULTED_4Q(cut+1:T+cut,:);
            Transfer_DEFAULTED_12Q=Transfer_DEFAULTED_12Q(cut+1:T+cut,:);
            bSIM_Defaultedm1=bSIM_Defaultedm1(cut+1:T+cut,:);
            
            bSIM_Defaultedm2=bSIM_Defaultedm2(cut+1:T+cut,:);
            bSIM_Defaultedm3=bSIM_Defaultedm3(cut+1:T+cut,:);
            bSIM_Defaultedm4=bSIM_Defaultedm4(cut+1:T+cut,:);
            bSIM_Defaultedm5=bSIM_Defaultedm5(cut+1:T+cut,:);
            
            BGSIM_Defaultedm1=BGSIM_Defaultedm1(cut+1:T+cut,:);
            BGSIM_Defaultedm2=BGSIM_Defaultedm2(cut+1:T+cut,:);
            BGSIM_Defaultedm3=BGSIM_Defaultedm3(cut+1:T+cut,:);
            BGSIM_Defaultedm4=BGSIM_Defaultedm4(cut+1:T+cut,:);
            BGSIM_Defaultedm5=BGSIM_Defaultedm5(cut+1:T+cut,:);

CA_SIM=CA_SIM(cut:T+cut-1,:);
SpreadSIM=SpreadSIM(cut:T+cut-1,:);
BadStandindSIM=BadStandindSIM(cut:T+cut-1,:);
priceSIM=priceSIM(cut+1:T+cut,:);   
zoneSIM=zoneSIM(cut+1:T+cut,:);   

VSIM=zeros(T,N);
%cSIM=zeros(T,N);
bindSIM=zeros(T,N);

TightSIM=zeros(T,N);
NUSSIM=zeros(T,N);
spreadSIM=zeros(T,N);
priceSIM=zeros(T,N);

for i=1:T
    for j=1:N
    bgrid_temp=squeeze(bgrid(BGiSIM(i,j),:)); 
    VSIM(i,j)=interp1(bgrid_temp,squeeze(V(BGiSIM(i,j),:,S_index(i))),bSIM(i,j),'linear','extrap');

 
   
    if DefSIM(i,j)==0
        spreadSIM(i,j)=(delta-delta*QSIM(i,j))/QSIM(i,j)-rbar;
        
    else
    end
    priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
    DebtLimSIM(i,j)=(KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))/qSIM(i,j);

   if abs(bpSIM(i,j)-(1/qSIM(i,j))*KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))<1e-3
            bindSIM(i,j)=1;
            TightSIM(i,j)=bpSIM(i,j)-(1/qSIM(i,j))*KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j));          
   end

    end
    
end
% Computing Solutions

CSIM=(omega*cSIM.^(-ita)+(1-omega)*yNSIM.^(-ita)).^(-1/ita);

meanY=mean(priceSIM.*yNSIM+ySIM);
Y_SIM=priceSIM.*yNSIM+ySIM;
% SUDDEN STOPS
CAo_SIM=CA_SIM./Y_SIM;  % Current Account
% Current Account and SS
SS=zeros(T,1);
SS_DE_th=(mean(CAo_SIM)-2*std(CAo_SIM));
SS(CAo_SIM<SS_DE_th)=1;


mean(SS);
BGSSm3_SIM=NaN(T,N);
BGSSm2_SIM=NaN(T,N);
BGSSm1_SIM=NaN(T,N);
BGSS_SIM=NaN(T,N);
BGSSp1_SIM=NaN(T,N);
BGSSp2_SIM=NaN(T,N);
BGSSp3_SIM=NaN(T,N);


bSSm3_SIM=NaN(T,N);
bSSm2_SIM=NaN(T,N);
bSSm1_SIM=NaN(T,N);
bSS_SIM=NaN(T,N);
bSSp1_SIM=NaN(T,N);
bSSp2_SIM=NaN(T,N);
bSSp3_SIM=NaN(T,N);



for i=1:T-6
    for j=1:N
        if SS(i+2,j)==1 
        BGSSm3_SIM(i,j)=BGSIM(i,j);
        BGSSm2_SIM(i,j)=BGSIM(i+1,j);
        BGSSm1_SIM(i,j)=BGSIM(i+2,j);
        BGSS_SIM(i,j)=BGSIM(i+3,j);
        BGSSp1_SIM(i,j)=BGSIM(i+4,j);
        BGSSp2_SIM(i,j)=BGSIM(i+5,j);
        BGSSp3_SIM(i,j)=BGSIM(i+6,j);

        bSSm3_SIM(i,j)=bSIM(i,j);
        bSSm2_SIM(i,j)=bSIM(i+1,j);
        bSSm1_SIM(i,j)=bSIM(i+2,j);
        bSS_SIM(i,j)=bSIM(i+3,j);
        bSSp1_SIM(i,j)=bSIM(i+4,j);
        bSSp2_SIM(i,j)=bSIM(i+5,j);
        bSSp3_SIM(i,j)=bSIM(i+6,j);
        end
    end
end






mean_default=mean(DefSIM);
cte_debtfactor=delta/(1-(1-delta)/R) ;
amean_PubDebt=cte_debtfactor*mean(BGpSIM)/mean(priceSIM.*yNSIM+ySIM);
amean_PvdDebt=mean(bpSIM)/mean(priceSIM.*yNSIM+ySIM);
% Private Crisis
BQSIM=(BGpSIM-(1-delta)*BGSIM).*QSIM;
avSpr=mean(spreadSIM); 
PvDebt_to_GDP=bpSIM./(priceSIM.*yNSIM+ySIM);
 

PubDebt_to_GDP=cte_debtfactor*BGpSIM./(priceSIM.*yNSIM+ySIM);
avPvD=mean(PvDebt_to_GDP);
avPub=mean(PubDebt_to_GDP);
NFA=PubDebt_to_GDP+PvDebt_to_GDP;

            BGSIM_Defaulted_DE=cte_debtfactor*BGSIM_Defaulted(~isnan(BGSIM_Defaulted))/meanY;
            bSIM_Defaulted_DE=bSIM_Defaulted(~isnan(BGSIM_Defaulted))/meanY;
bSIM_Defaultedm1_DE=bSIM_Defaultedm1(~isnan(BGSIM_Defaulted))/meanY;
bSIM_Defaultedm2_DE=bSIM_Defaultedm2(~isnan(BGSIM_Defaulted))/meanY;
bSIM_Defaultedm3_DE=bSIM_Defaultedm3(~isnan(BGSIM_Defaulted))/meanY;
bSIM_Defaultedm4_DE=bSIM_Defaultedm4(~isnan(BGSIM_Defaulted))/meanY;
bSIM_Defaultedm5_DE=bSIM_Defaultedm5(~isnan(BGSIM_Defaulted))/meanY;


BGSIM_Defaultedm1_DE=cte_debtfactor*BGSIM_Defaultedm1(~isnan(BGSIM_Defaulted))/meanY;
BGSIM_Defaultedm2_DE=cte_debtfactor*BGSIM_Defaultedm2(~isnan(BGSIM_Defaulted))/meanY;
BGSIM_Defaultedm3_DE=cte_debtfactor*BGSIM_Defaultedm3(~isnan(BGSIM_Defaulted))/meanY;
BGSIM_Defaultedm4_DE=cte_debtfactor*BGSIM_Defaultedm4(~isnan(BGSIM_Defaulted))/meanY;
BGSIM_Defaultedm5_DE=cte_debtfactor*BGSIM_Defaultedm5(~isnan(BGSIM_Defaulted))/meanY;

           
            DistSIM_DEFAULTED_4Q_DE=DistSIM_DEFAULTED_4Q(~isnan(BGSIM_Defaulted))/meanY;
             DistSIM_DEFAULTED_12Q_DE=DistSIM_DEFAULTED_12Q(~isnan(BGSIM_Defaulted))/meanY;
            Transfer_DEFAULTED_4Q_DE=Transfer_DEFAULTED_4Q(~isnan(BGSIM_Defaulted))/meanY;
            Transfer_DEFAULTED_12Q_DE=Transfer_DEFAULTED_12Q(~isnan(BGSIM_Defaulted))/meanY;


BGSSm3_DE=cte_debtfactor*BGSSm3_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSSm2_DE=cte_debtfactor*BGSSm2_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSSm1_DE=cte_debtfactor*BGSSm1_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSS_DE=cte_debtfactor*BGSS_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSSp1_DE=cte_debtfactor*BGSSp1_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSSp2_DE=cte_debtfactor*BGSSp2_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSSp3_DE=cte_debtfactor*BGSSp3_SIM(~isnan(BGSSm3_SIM))/meanY;


bSSm3_DE=bSSm3_SIM(~isnan(bSSm3_SIM))/meanY;
bSSm2_DE=bSSm2_SIM(~isnan(bSSm3_SIM))/meanY;
bSSm1_DE=bSSm1_SIM(~isnan(bSSm3_SIM))/meanY;
bSS_DE=bSS_SIM(~isnan(bSSm3_SIM))/meanY;
bSSp1_DE=bSSp1_SIM(~isnan(bSSm3_SIM))/meanY;
bSSp2_DE=bSSp2_SIM(~isnan(bSSm3_SIM))/meanY;
bSSp3_DE=bSSp3_SIM(~isnan(bSSm3_SIM))/meanY;

BG_goodtimesSIM=BGSIM(CAo_SIM>=0);

fprintf('The long-run probability of crisis in DE is: ');
fprintf('%5.4f \n ',  mean(SS)*100);


fprintf('The long-run probability of being in the shrinking default zone in DE is: ');
fprintf('%5.4f \n ',  mean(zoneSIM)*100);


DefSP_DE=DefSP_SIM(~isnan(DefSP_SIM));

fprintf('Total number of defaults in DE is: ');
fprintf('%5.4f \n ',  size(DefSP_DE));

fprintf('Total number of defaults in DE that planner would have avoided is: ');
fprintf('%5.4f \n ',  sum(DefSP_DE));


meanYDE=meanY;
PubDebt_to_GDP_DE=PubDebt_to_GDP;
PvDebt_to_GDP_DE=PvDebt_to_GDP;
Def_funcDE=Def_func;

